Multiconfigurational time-dependent Hartree method for describing particle loss due to absorbing 

boundary conditions 



Simen Kvaal 12 E] 

1 Matematisches lnstitut, Universitdt Tubingen, Auf der Morgenstelle 10, D-72076 Tubingen, Germany 
2 Centre of Mathematics for Applications, University of Oslo, N-0316 Oslo, Norway 

Absorbing boundary conditions in the form of a complex absorbing potential are routinely introduced in 
the Schrodinger equation to limit the computational domain or to study reactive scattering events using the 
multi-configurational time-dependent Hartree method (MCTDH). However, it is known that a pure wave- 
function description does not allow the modeling and propagation of the remnants of a system of which some 
parts are removed by the absorbing boundary. It was recently shown [S. Selst0 and S. Kvaal, J. Phys. B: 
At. Mol. Opt. Phys. 43 (2010), 065004] that a master equation of Lindblad form was necessary for such a de- 
scription. We formulate a multiconfigurational time-dependent Hartree method for this master equation, usable 
for any quantum system composed of any mixture of species. The formulation is a strict generalization of 
pure-state propagation using standard MCTDH. We demonstrate the formulation with a numerical experiment. 



I. INTRODUCTION 

Today, the de facto standard approach in ab initio quantum- 
mechanical many-particle propagation is the multiconfigura- 
tional time-dependent Hartree method (MCTDH) and its vari- 
ations QJ0. Already forN = 2 electrons in three dimensions, 
the full six-dimensional time-dependent Schrodinger equation 
is very hard to solve and can only be handled on supercom- 
puters. With MCTDH for identical particles the exponential 
scaling of the Hilbert space dimension with respect to N is 
"postponed" to higher number of particles, and the N = 2 
propagation can be done on a single desktop computer. Cur- 
rent implementations can handle Af < 8 electrons in cylindri- 
cal geometries reliably |8]|9). For bosons, the Pauli exclusion 
principle is absent, and MCTDH can treat hundreds ifTUl and 
even thousands of particles ifTTl in one-dimensional geome- 
tries, and recent multi-layer MCTDH techniques allow distin- 
guishable dimensions in the thousands with relative ease |7j. 
showing great promise of extending the domain of application 
of MCTDH methods. 

MCTDH is derived using the time-dependent variational 
principle lfl2l [PJl . As such, it is energy conserving, unitary 
and quasi-optimal in the sense that the growth of the error in 
the 2-norm is locally minimized. 

Ab initio dynamical problems in quantum mechanics are 
formulated on an infinite domain which must be truncated for 
numerical calculations. The numerical reflections implied by 
the truncations are usually dealt with using absorbing bound- 
ary conditions of some sort. The most common approach is 
to introduce a complex absorbing potential (CAP) in a region 
around the truncated domain fl4l . That is, the Hamiltonian H 
is mapped to H — iF, where T > is a local one-body potential 
vanishing on the domain of interest, and only taking nonzero 
values outside the domain. This approach is also used in or- 
der to calculate properties like reaction and ionization prob- 
abilities, and CAPs are routinely implemented in MCTDH 

codes he nana. 
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Other absorbing operators are also common, such as the 
so-called transformative CAP (TCAP) lfl8l . which is more or 
less equivalent with the non-local CAP obtained using smooth 
exterior scaling [19] or perfectly matched layers [20 1. While 
exact and space-local absorbing boundary conditions may be 
formulated [21], they are in general non-local in time, and 
therefore impractical. In this work, we focus solely on a local 
CAP for simplicity, but any absorbing operator can be used. 

Given a system of N particles, the wave-function is nor- 
malized to the probability of finding all particles within the 
computational domain. With a CAP, evolves according to 
the non-Hermitian Schrodinger equation 

i^sf N = {H-sysf N . (i) 

An elementary calculation gives ajH^jvH 2 = -2( l P w |r| l P w ) 
for the probability derivative. Consequently, if the wave- 
function overlaps the CAP, the whole wave-function decays 
and eventually vanishes; it does not approach a wave-function 
with a different number of particles. In other words, even with 
an absorbing boundary, one is stuck with an Af -particle de- 
scription. Information like ionization probabilities and reac- 
tion rates may be obtained from evolving Eq. ([T]) alone, but if 
the remainder of the system is desired, i.e., a description of 
the N — I, N — 2 particle systems, etc, one is at loss. 

In a recent article it was argued that the solution is a den- 
sity operator approach ll22l because the loss of particles is 
an irreversible process; H — iF is a non-Hermitian operator 
implying a preferred direction of time. The necessity of the 
quantum dynamical semigroup describing the evolution to be 
trace-preserving, Markovian and completely positive implies 
the applicability of the famous theorems due to Lindblad and 
Gorini and coworkers [23 -25], giving a master equation on 
Lindblad form. The resulting equation is 

j ( Pn = -i[H,Pn] - {r,p„} +2 J r(x)^(jc)p„ +1 ^(x) t dx, 

(2) 

where p„ is the density operator for the n-particle subsystem, 
< n < N. The integral is over all discrete and continous 
degrees of freedom. It should be noted that the non-Hermitian 



2 



Schrodinger equation ([T]) is equivalent to the von Neumann 
equation 



dt 



Pn 



-i[H,p n ]-{T,p n } 



for each n-particle system. The original formulation of the N- 
particle problem is not changed, but extended to yield n < N 
particle systems as a by-product. Importantly, nowhere is an 
ad hoc hypothesis introduced. 

In this article, we formulate a MCTDH method for Eq. (|2]i 
for identical particles and mixtures. It is based on a so-called 
type II density operator manifold [3, 26 1 and is fully varia- 
tional. It turns out to be exactly trace-preserving and a strict 
generalization of MCTDH for pure states evolving according 
to the non-Hermitian Schrodinger equation ([TJ. The method 
is first formulated for identical fermions and bosons, and then 
extended to arbitrary mixtures of species. The derivation for 
identical particles is done in Sec. [Ill] while the generalization 
is done in Sec. [V] after a numerical experiment on a system 
of identical fermions is presented in Sec. IV Not surprisingly, 



the resulting formulation is a direct generalization of the pure 
state MCTDH treatment of mixtures [27|. 

To prepare for the density operator MCTDH method (here- 
after called p-MCTDH as opposed to ^-MCTDH for pure 
states), we give a brief derivation of Eq. Q in Sec. [TT] that 
underlines the inevitability and uniqueness of the master equa- 
tion on Lindblad form. We also show how the probability in- 
terpretation of ll'IVl) 2 , which is not purely quantum mechan- 
ical, can be interpreted in terms of measurements performed 
continuously in time. 



II. LINDBLAD EQUATION FOR SYSTEMS WITH A CAP 
A. Classical probabilities from a CAP 

The evolution of under the non-Hermitian Schrodinger 
equation ([T| is irreversible. It introduces a preferred direction 
in time, which is easily seen from the fact that the evolution 
cannot be reversed: eventually W^nW > 1, and the backward 
propagation may be non-existent, even mathematically. 

We now give an interpretation of the square norm ||1V|| 2 
in terms of measurements performed continuously in time, 
thereby exposing the irreversibility. To this end, consider first 
a closed single-particle system, described by the Hamiltonian 
H — without a CAP — in the Hilbert space M\ . Quantum me- 
chanically, the squared norm H'Pill 2 of the wave-function is 
the probability of finding the particle somewhere in configura- 
tion space. Equivalently, it is the probability of obtaining the 
value 1 upon measurement of trivial observable / (the identity 
operator). 

Adding a CAP — ir to the Hamiltonian and keeping the 
probability interpretation of ||M'i|| 2 has non-trivial implica- 
tions. The description is obviously no longer in accordance 
with the basic postulates of quantum mechanics, since 



that is to say, total probability is not conserved. 

Suppose we perform a single measurement on the observ- 
able P given by 



|jc) (xjdjc, 



i.e., of the projection operator P onto £2, the truncated compu- 
tational domain. It has two eigenvalues, and 1, correspond- 
ing to finding the particle outside or inside £2, respectively. 
We will obtain the answer with probability 1 — H^iHq = 
1 - (*Pi l-Pl^i) and the answer 1 with probability H 1 ?! Af- 
ter measurement, the wave-function collapses onto (I — P^i 
in the former case, and P*¥\ in the latter. Note that the wave- 
function collapse is irreversible, as the original wave-function 
cannot be reconstructed after the event. 

Suppose we perform many experiments at short time inter- 
vals t = ht, n = 0, 1,2, . , .. Each experiment yields certain in- 
formation about the system after the measurement: we know 
with certainty if the particle is in £2 or not. After n experi- 
ments, all giving 1 as answer, the wave function is (to first 
order in r) 

S?i(irr) = [Pe- irH ]' !l Pi(0), 
and the probability that n positive answers have been given is 

ll^iMH 2 , 

again to first order. This probability is classical in the sense 
that it is a probability description of the history of our macro- 
scopic measurement device, or of its print-out on a sheet of 
paper if one prefers. 

If we allow the approximation 

P ss exp(--zT), 

which dampens 'J'i strongly outside £2, we get 



Pe 



AtH , 



-tF-itH _ 



Ar(H-ir) , ft 2 ) 



dr 



-2<^i|r|^i)<0, 



In the limit of small r, we see that the solution to the 
Schrodinger equation with a CAP is obtained. ||Wi(f)|| 2 is 
then the probability of finding the measurement apparatus in 
the state "the particle has not yet been found outside £2" at 
time t. 

Our discussion is immediately suggestive of interpreting 
the CAP as an actual model for some external detecting de- 
vice, which was already pointed out in the first paper on CAPs 
in quantum mechanics fl4l . Although it is likely that Eq. |2]) 
can be derived in such a way, it is not relevant here, as the 
Schrodinger equation with a CAP is our starting point. The 
above discussion is only intended as a means for understand- 
ing the irreversibility of the non-Hermitian dynamics. What- 
ever interpretation we use, we see that the single-particle sys- 
tem undergoes a transition to a zero-particle state (i.e., zero 
particles in £2) in an irreversible way, and that the probability 
of having n = 1 or n = particles is not quantum mechanical. 

A complete specification of the quantum state of the n < 
1 -particle system is a density operator in the Hilbert space 
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M\ ®J%b, where Jf„ is the Hilbert space of n particles. In 
this case, the density operator is, on block form, 



WFil 

o (i-ofi!^)) 



Note that M{) is one-dimensional, and the lower right block is 
just a non-negative real number, the probability of zero parti- 
cles in il. 

For a rigorous treatment of the above limiting process, see 
Chapter 7.4 in Ref. EH). 

Generalizing the discussion to N particles is completely 
analogous to our discussion, approximating P s» exp(— ir) in 
the same way. ||\PjvJ| 2 is then the probability of the measure- 
ment apparatus showing that all particles are inside £1 



B. A'-fermion systems with a CAP 

Suppose a complete set of fermionic or bosonic creation 
operators {cj}, j = 1,2,..., are given. An important formal 

property of c T j is that it defines a map from to the 
Hilbert spaces of n and n + 1 fermions, respectively, and that 
they fulfill the (anti-)commutator relation 



{c],q-}± = c]q-±qc] = Sj k 



(3) 



with the plus sign for fermions, and the minus sign for bosons. 
For us, it is natural to work within Fock space Jf?, defined 

by 

CO 

«=0 

containing all possible states with any number of particles. As 
operators in Jf, the ci are orthogonal in the Hilbert-Schmidt 
inner product, 



and traceless, 



[(cj,c k ))=tr(c)c k )=0, j^k, 



tr(c ; -)=0. 



Given an n-particle Hamiltonian on first-quantized form, 
viz, 

H n = th(i)+"t t 

i=l i=l j=i+l 

(where the indices i and j in the operators indicate which par- 
ticles' degrees of freedom they act on,) we may write it com- 
pactly on n-independent form using the expression 



jk jklm 



^c\c m ci, (4) 



where fin and um m are the usual one- and two-particle inte- 
grals, respectively. 



Adding a cap F — YJ!=i^(i) amounts to modifying the 
single-particle coefficients since 



(5) 



It is important to note, that given the set {cj}, the second- 
quantized expressions of H and F are unique, and vice versa. 



C. The Lindblad equation 

We now discuss the master equation introduced in 
Ref. Il22l . for a system of identical particles. We underline that 
the master equation follows uniquely from the probability in- 
terpretation of H^nll 2 for an n-particle wave-function, and the 
requirement of a Markovian, trace-preserving and completely 
positive evolution. We also address some mathematical points 
omitted in Ref. El . 

Describing a system with a variable number of particles, we 
work with states in Fock space. Our starting point, the non- 
Hermitian Schrodinger equation with a given CAP, then reads 



i— l P= (H-iFW. 

At 1 ' 



(6) 



With an initial condition with exactly particles, i.e. ^(O) = 
G '^Ni this equation is equivalent to Eq. ([TJ, which we 
repeat here for convenience: 



i-v N = (H-nyv N . 



(7) 



This description is valid for all particle numbers N. 

As Eq. (|6|l is irreversible, it is not possible to find some new 
Hamiltonian H' in Fock space that generates a unitary evo- 
lution which describes a decreasing number of particles and 
also reproduces Eq. (17). Instead, one must turn to a density 
operator description. For a density operator p € TC(Jtff), the 
trace-class operators in Fock space l25ll28l . Eq. (|6]l is equiva- 
lent to the non-trace conserving von Neumann-equation 



-i[H,p]-{F,p}, 



(8) 



which is verified by a simple computation. The density oper- 
ator p has a natural block structure. Letting P„ be the orthog- 
onal projector onto M^, we have the resolution of the identity 

n=0 

Applying this to either side of p gives the decomposition 

P = ^ Pn.ni'. Pnjn — PnpPm- 
n=0m=0 

This block structure is depicted in Fig. [T] Projecting Eq. ([8]) 
from the left and right with P n , we obtain 



dr 



pn,n = -i[H,p n ,n] ~ {T,Pn,n}, 



(9) 
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FIG. 1: Illustration of the block structure of the density operator 
in Fock space with at most N identical particles. In general, each 
block is infinite dimensional (except for po which is 1 x 1). The 
diagonal blocks are emphasized, as for pure state initial conditions 
P = I^n) O^ivl me density operator turns out to be block diagonal. 



which is equivalent to Eq. |7]). 

We wish obtain a master equation for p that does preserve 
total probability. At all times, there is a probability tr(p„) of 
having n particles (we define p„ = p n ^), and these should add 
up to tr(p) = L„tr(p„) = 1. Since the von-Neumann equation 
is local in time, we require the master equation to be Marko- 
vian. Our master equation must reproduce (|9|l for n = N, 
where N is the number of particles initially present. Of course, 
foin<N the master equation will end up with some source 
terms. This must hold for all N, i.e., our master equation 
should not depend on N; only the initial condition depends 
on N. 

We also require the master equation to describe a physi- 
cal system which in turn may interact with other quantum 
systems. A minimal requirement for such master equations 
is that its quantum semigroup is completely positive. Com- 
plete positivity roughly states that if our system is described 
together with a different system with Hilbert space "V , so that 
the combined state is <x g TC(^€5 Y), then the semigroup 
of cr preserves the self-adjoint positive-semidefiniteness of a. 
Surprisingly, requiring this property on the flow of p alone is 
not enough E4ll25l . 

Let us for the moment assume that Fock space has finite 
dimension, i.e., it is defined by finite number L of creation op- 
erators. The theorem by Gorini and coworkers is now appli- 
cable |24|: Any trace-preserving, Markovian and completely 
positive quantum dynamical semigroup has a master equation 
on the form 



-p = -i[H,p]+S)(p), 
where the dissipative terms have the generic form 



(10) 



@[p) = Y.dapQLppLl - {4W). (H) 

The operators L a form a traceless, orthogonal set with the 
Hilbert-Schmidt inner product ((p,cr)) = tr(p'<x) (the operators 



are then linearly independent in TC(J^)), and [d a p] is a Her- 
mitian positive semi-definite matrix. Importantly, the quan- 
tum dynamical semigroup generated by Eq. ( [T0| is uniquely 
given by H and 3 and vice versa. 

The L a play the role of a basis in TC(Jtff), which is of di- 
mension dim(Jf ) 2 . It is of course not unique, but given a 
choice of basis, the coefficients [d a p] are unique. Notice that 



the trace of the anti-commutator term in Eq. ( 1 1 1 is equal to 
but of opposite magnitude compared to the term where p is 
sandwiched between Lp and L\. 

The von-Neumann equation is almost of the form 



10 1; it 



lacks the "sandwiched terms" responsible for a compensation 
of the trace decrease due to the anti-commutator, i.e., 



p = -i[H,p\-Yjjk{c)c k ,p\. 

jk 



As the operators {c,} are indeed traceless and orthogonal, it is 
immediately clear that the correct master equation is obtained 
by simply adding the sandwiched terms, i.e., 



p = -i[H,p] -Y^jk{c)ck,p}+2Y^jkCkPc). (12) 



At 



If the dissipative terms are not chosen on exactly this form, 
Eq. Q will not be reproduced for all N. 

Indeed, consider an initial condition on the form 



P (o) = \y N )(v N \. 



(13) 



Projecting Eq. ( p"2} from the left and right with P n and P m , 
respectively, we obtain the differential equation obeyed by the 
block p n<m , 



dt 



Pn,m = -i[H,Pn,m\ - {T,Pn,m} + ^Yj~jkCkPn 



+ l,m+lC 



(14) 

We see that the off-diagonal blocks p n m with n ^ m are iden- 
tically zero with the initial condition ( p"3j ). The compensat- 
ing sandwich terms are seen to be responsible for transporting 
probability from the -particle system into the — 1 parti- 
cle system, and so on, downwards along the diagonal blocks 
in Fig. [T] Moreover, the evolution of pAyv is equivalent to 
Eq. (|7ji. Modifying the coefficients Y ^ or adding more lin- 
early independent operators to the set {c,-} will void this for 
some N. 

The theorem by Gorini and coworkers is only valid in 
a finite dimensional Hilbert space. For infinite dimen- 
sional spaces, Lindblad discovered a generalization for norm- 
continuous semigroups ||23l , for which H and Y are necessar- 
ily bounded operators. This is rarely the case, but a Lindblad- 
type theorem for the general unbounded cases is simply not 
known |25|. On the other hand, all known examples of 
completely positive semigroups in infinite dimensional spaces 
have generators on the generic form 
bounded H or [d a p\ . 

It now seems reasonable to remove the restriction of a 
finite-dimensional Fock space, which does not change the for- 



1 1 1, with possibly un- 
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mal appearance of the Lindblad equation ( [12) , Using the rela- 



tion 



j 

where {<Pj(x}} is the orthonormal single-particle basis func- 
tions associated with {cj}, we arrive at the master equation 

j t P = -i[H,p] - {Up} +2 J r(x)^(x)p^r(x) t dx, (15) 
which on block form becomes 

d p„ = -i[H,p n ) - {r,p„} +2 [ r(x)^(x)p n+i i/f(xy dx. 



dt 

(16) 

where we have assumed an -particle initial condition p(0) = 
Pn- 



We stress that the Lindblad equation (16i followed only 
from the probability interpretation of H^jvjp anc ' fr° m the 
requirement that the master equation generates a Markovian, 
trace-preserving and completely positive semigroup. 



III. MCTDH FORMULATION FOR IDENTICAL 
PARTICLES 

A. The p-MCTDH manifold 

We now derive the p-MCTDH approximation to the Lind- 



blad equation ( 15 1. This method is in a sense a combination 
of "P-MCTDH for identical particles Rti6ll27ll. using second 
quantization, and p-MCTDH |3 26|. Our method is necessar- 
ily formulated in Fock space and describes a variable number 
of particles. To the best of our knowledge, p-MCTDH has 
not been formulated using second quantization in the litera- 
ture. Moreover, the global variational principle employed in 
"P-MCTDH (the action integral point of view) is not appli- 
cable for dissipative systems, since 2#(p) is not self-adjoint. 
The second-quantization techniques used here are therefore 
somewhat different from the pure-state approach. As the par- 
ticular properties of the equations of motion for the present 
problem are more easily exposed in a thorough investigation, 
we choose to do a detailed derivation. 

We repeat the Lindblad equation ( p~5] > for convenience: 



dt 



P = J?{p) 
= -i[H,p] + 9(p) 
= -i[H,p]-{T,p}+2 



T(x)4f(x)p^(x) f dx. (17) 



Here, p g TC(Jf) is a Fock space density operator, which (in 
the exact equation) is block diagonal. ifr(x) destroys a particle 
at the configuration-space point x € X. Typically X = K 3 x 
{t,i} for a spin-1/2 fermion, i.e., x = (r,m), although our 
derivations are completely independent of X. 

For all t > 0, p(f) is approximated by an element in a 
manifold j$ C TC(J^), with inner product inherited from 



TC(^#). (Strictly speaking, the inner product is inherited 
from the Hilbert space of Hilbert-Schmidt operators, in which 
TC(J%?) is dense. Not all Hilbert-Schmidt operators are trace- 
class.) An approximate variational differential equation on ^# 
is sought. The time-dependent variational principle ll2l [T2l[T3ll 
chooses the time derivative p = dp/dt to minimize the local er- 
ror in the norm induced by ((■, •)) as follows: 



((6p,p - &(p))) = 0, V(5p G T P J{, 



(18) 



where T p jj( is the tangent space at p, i.e., the space of all 
possible time derivatives of p. Consequently, the right-hand- 
side Jz?(p) of Eq. ( 17 1 is projected orthogonally onto T p ^t, 
and we have 



argmin | 



■JZp\\ 



We choose ^# as the so-called "type II" density operator 
manifold (3j [17] [26], albeit with a slight generalization as 
we consider Fock space instead of a fixed number of parti- 
cles. There is an alternate way of defining a variational mani- 
fold, the "type I" manifold, but it does not reduce to the usual 
MCTDH in the case of a pure state. 

The manifold is defined as follows. Given a finite set <p 
of L single-particle functions (SPFs) ipj E J%i, 1 < ./' < L and 
their corresponding creation operators cj, 



cpj(x)iff (x) f dx, 



we consider the subspace % of spanned by all possible 
linearly independent n-body functions built using products of 



The notation J[n] means an ordered tuple (juj2>'" >jn) of n 
single-particle indices, i.e., jfi < j% < ■ ■■ < j n . For fermions 
&j[ n ] is a Slater determinant and ji < }2 < ■ ■ ■ < j n , while for 
bosons it is a permanent. We then consider the subspace Y of 
Fock space, 

N 

r = (&% 

n=0 

spanned by all the <t> j, where J means a ordered tuple with n 
indeterminate. As we describe at most particles, we trun- 
cate at n < N in order to ensure a finite-dimensional space. 
(For bosons any number of particles may occupy each ipj, cre- 
ating an infinite dimensional space if we do not truncate the 
sum, even with L finite.) 

Each p G is now defined as an arbitrary linear operator 
in Y, viz, 



--^j)Bjk(^k\ 

JK 



We see that p is a matrix with respect to a time-dependent 
orthonormal basis. We denote by B the matrix formed by the 
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Bjk, i.e., the Galerkin matrix. To sum up, p is parameterized 
in terms of an arbitrary matrix B with respect to the basis gen- 
erated by an arbitrary set of L SPFs tp. 

The set tp may formally be extended to a complete ONB 
tp for M\, a (usually infinite) set of functions tp s , s> L, such 
that the second-quantized Hamiltonian is given by Eq. Q, and 
the CAP by Eq. |5]), but where the expansion coefficients in 
general will depend on the particular value for tp and tp. This 
will be of use later on. 



B. Parametric redundancy and tangent space 

For a given p £ the parameters tp and B are not unique. 
Since "V is determined only by the subspace spanned by tp, not 
the individual tpj, any unitary change 



k 



G = U(L) 



where U(L) is the unitary group of L x L matrices, yields the 
same space "f, and therefore the same operators p can be pa- 
rameterized. Under the group element G, the basis functions 
transform as 



4> 



Jn] 



fcl k„ 



where the sum Ljf [„] is over a ^ multi-indices of length n, and 
not only ordered ones. Defining the transformation of B by 



Bjk 



J'K' 



we see that 



P 



■P- 



Moreover, B jk are all independent parameters, showing that 
U (L) is in fact the largest group of transformations leaving p 
invariant. 

The non-uniqueness of tp and B implies that given a deriva- 
tive (tangent vector) p £ T p ^0, the derivatives tp and B are 
not unique. Suppose p(f) £ j% is a given smooth path. There 
exists tpo(t) and Bq(?) such that p(t) = p(tpo(t),Bo(t)). By 
the considerations above, any other possible parameter path 
is on the form (tp(t),B(t)) = (ipo(t)G(t),@(t) f Bo&(t)), where 
<po(t)G(t) stands for the transformation 

foj(t) — ► <pj(t) = Y,<Po,k{t)G k j(t). 

k 

(We consider tp a "matrix" whose columns are tpj.) All deriva- 
tives of tpo are of the form 

£o = W7 +X> (fojWk) = 0. 

The functions Xj are all independent. Since GG^ = II, we find 
that G = gG, with — g^ = g £ u(L), the Lie algebra of the Lie 
group U(L). The transformed tp = tpoG then has the derivative 



dr 



tp G= [tp (i] + g) +x\G, 



and it is seen that if we choose g = —rj, then tp is in fact unique, 
since G = r\G uniquely specifies G(t). This is equivalent to the 
condition 

(<Pj\tp k )=0, Vj,k. 

In this way, there is a there is a one-to-one map between 
triples (tp,B,g) and p, with g £ u(L). The element g is then 
called a gauge choice, and the gauge choice induces a unique 
parameterization (tp(t),B(t)) of p(f). This kind of differential 
geometrical structure is called a principal bundle |29|, and is 
familiar in quantum field theory — but it arises in a completely 
different way! 

It is easily verified that 

where D is the operator 



j=l 

and where ct is as the operator that creates the single-particle 
function tpj. We observe that for any two single-particle func- 
tions u and v, not necessarily normalized, the relation 

{c{u),c\v)} ± EE C(u)c\v)±c\v)c(u) = (u\ V ) 

is obtained by expanding each operator in the field creation 
operators. 

We are now ready to consider an arbitrary time derivative 
of an element p(t) £ 



p = £ \<i>j)B lK (<P K \ + \<Pj)B lK (<P K \ + \<S>j)B JJC (<Z> K \ 

J,K 

= Dp+pD f + ^\<Pj)B lK (<P K \. (19) 



J.K 



In order to perform the projections in Eq. < pT8j >, we must 
identify all linearly independent tangent vectors, i.e., all inde- 
pendent admissible infinitesimal variations of p. This amounts 
to varying the B jk independently, and the tpj independently, 
but according to a specific choice of gauge. For simplicity, we 
consider the choice (tpj\tpk) = 0, which generates the simplest 
equations and are also the most common in MCTDH theory. 

From this it follows that the admissible time derivatives of 
tpj are arbitrary functions # = Q&, where 



Moreover, 



Q=i-YM){vk\ 

k 



(*/|*jc) = <*/|Z>|**)=0, VJ,K. 



Now, the independent variations of p can be divided into 
two groups: For each pair J, K, the matrix element Bjk can 
be changed, giving a tangent vector 



«sp = |*/><**l 



(20) 



7 



An arbitrary change & — Q§ in tpj consistent with the gauge 
choice gives 



5p = C* (p)c jp + pcjcfo) . 



(21) 



Inserting these two expressions into the variational principle 
will yield a complete set of differential equations for B in the 
first case, and <p in the latter. 



C. Equations of motion 

We use the notation to indicate a sum over the complete 
set of SPFs. We let H be the Galerkin matrix of H, i.e., 

H JK = (<Pj\H\<P K ), 

and analogously define G to be the Galerkin matrix of F. We 
let Cj be the Galerkin matrix of cy, which in fact is independent 

of <p: 



(cj)jk = (*j\cj\<S> K ) = (VjWjWk 



(22) 



The primed quantities corresponds to any other choice of 
SPFs. The independence follows from the (anti-)commutator 
Q which only depends on orthonormality. 

The Galerkin matrices of H and T can be expressed as 



H = Y,K<P)jk c ) c k + \ £ u(<p)jki m c]c\c m ci (23) 

(24) 



' jkhn 



G = £r(^cjc, 

jk 



The expansion coefficients are dependent on <p at time f, but 
the creation- and annihilation matrices are not. The Galerkin 
matrices are naturally expressed using some fixed, abstract 



basis due to Eq. (22i. Existing methodology for computing 



matrix-vector and matrix-matrix products can be re-used by 
referring to this basis. 

To derive the equations of motion, we begin by inserting dp 



from Eq. (|20[> into the variational principle (18 1. For the term 



((<5p,p)) and the term containing [H,p], and for the sandwich 
term, we obtain, respectively, 



|*/)<**|, £ bj>K>\*f)(*K'\+Dp+pDty =tr 



J'K' 



\®k)(<5>j\ £ \<S>j,)Bj, K ,{<S> Kl \+Dp+pD 



l^^l.-i^pj^-itr 



1^)^/1 ^£ \®j')Bj>k'{®k>\- £ \<S>j>)Bj,v{<i> K ,\H 



'\<Pj)(<P K \,YF jk c k pc)))=tr 



jk 



\*k)(*j\ £r # Q £ \<t> r )Bj, K ,(<i> K ,\c} 

\ jk J'K' , 



= Y / Tjk(CkBc 1 j)j K . 

jk 



Bjk (25) 
l-iH,B])« (26) 
(27) 



In Eq. (25 1 the terms containing D vanish since (pj\(pk 
{cj,c\}± — 0. A similar calculation as Eq. (26 1 yields 



I*/)(*4{T.P}) = ({G,B})«. 



Assembling Eqs. ( |25} through ( |28| , we get the equation of 
motion for B: 

B = -i[H,B]-{G,B}+2£r,*c jfc Bct. 

jk 



We now make the observation that B^ = B, showing that p^ = 
p is preserved during evolution, which we may use when we 



turn to the projection onto the tangent vector in Eq. (21 1. Let 
(28) F be an arbitrary operator, and calculate 



((6p,Fp+pF f )) = tr |[c()?) t c j p+pc ; t c(^)](Fp+pF t )} = tr c{&) t c j pFp + c{§) f c j p 2 F 1i +pc)c(§)Fp+pc]c(V)pF J < 
= 2Retr[c]c(i?) J Fp 2 ], 

I 



since c(§)p = 0, as c(§) annihilates a function orthogonal to Setting F = D = Y,k^\ c k we obtain 
all the <pj. 

((6p,Dp+prf)) = £2Retr[c) C ( 1 ?)c A t c it p 2 ] 

k 

= £2Re(^|^)tr(4 Q p 2 ), 

k 
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and we note in passing that ((dp, \^j)Bjk(^k\ 
since c(ff)\&j) = 0. 

The "sandwich" term in the master equation also gives zero 
contribution, since 

tr [c(^c jP c lP ci] = tr\pcic{$) f c jP } = 0, 

since c(#)c m \$j) = 0. 

In these calculations § = Q§ was arbitrary. Choosing — i§ 
instead turns "Re" into "Im", so we may drop taking the real 
part. Assembling this, we get the equation 



0, again Upon insertion in Eq. ( 29 1, we find 



vY,m k Mc)c k p 2 )=tx[c)c^)Fp 2 



(29) 



tx[c)c(#)Fp 2 ] = ^.fk,^[c)c{M c iP 2 \ 

kl 

= ZZ(^){<Pk\f\<PlMc}ciP 2 ) 



k I 

= J>|/|w)tr(cJqp^ 
/ 



In the last calculation, we used Y,jWj) (<Pj\ = I< 
Second, suppose F is a two-particle operator, viz, 



with F = H — ir, and the equation must hold for all & and all 
J- 



We now compute the right-hand-side of Eq. ( 29 1 for arbi- 
trary single-particle and two-particle operators F. Suppose at 
first 

F = Y/j kC ) Ck ' fa = (vjlfM- 

jk 



F = 9 H fjklmc)c k C m C h f jklm = {<Pj<Pk\f\<Pl<Pn 



' jklm 



where we assume fj k \ m — fkjmh The inner product is on 3%{ eg) 
i.e., the brackets are not antisymmetrized. Now the right- 



hand-side of Eq. ( 29 1 becomes 



tr[c'jC{&)Fp 2 } = ~ Y,fklrm^[c)c{ff)c\c}c n c m p 2 ] = £ (<Pk<Pi\f\<Pm<P„) tr{c]c] c„c m p 2 ) 



' klmn 



k bnn 



Imn 



r 



In the first step, we used the symmetry of fkimn ar, d the (anti- 
)commutator relation for the creation operators operators, the 
result being the same regardless of particle statistics. 

Three-particle operators, or even higher, are computed in 
similar fashion. For a three-body operator 

F = ^^^fjklpqrCjclcJcrCqCp 
J ■ jklpqr 

we obtain the right-hand side 
lr[c)c{u)Fp 2 ] = £ {fty m \f\y p ip q tp,)tr(c)clc}c r c q c p p 2 ) 

klpqr 

In all cases, the combinatorial factor cancels due to symmetry 
properties. 
The matrix 



S jk = tr(c)c k p 2 )=tr{c)c k B 2 



(30) 



defines the p-MCTDH analogue of the reduced one-body 
density matrix entering at the same location in standard \P- 
MCTDH theory. Similarly, the analogue of the reduced two- 
body density matrix is defined by 



S Tk'im = K(c)clc m c,p z ) = tr(c}cT Cm C/B / 



(31) 



and so on. 



We may now assemble the various one- and two-body con- 
tributions to the SPF equation of motion: 

k k klm 

which holds for all # — Q§. Since ipj = Qtpj, we arrive at the 
final single-particle equations of motion: 

i^foS /* = Qjl( h ~ iF )fkS j k + Q( ■ <p k \u\<pi<p m )<piS f^ lm 

k k klm 

= fi£(A-iI>*S j k + ^QUkm<PiS%, 

k klm 

where 

( ■ <Pk\u\<pi<p m ) = J (pk(y)u(x,y)ipi(x)tp(y) dy, 

and where the mean-field potentials U km are defined by 

U km (x) = J (Pk(y)u(x,y)<p(y) dy. (32) 

Assuming u(x,y) to be a local potential, U km (x) is also a local 
one-body function. 
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D. Discussion 

Let us sum up the equations of motion for the density oper- 
ator p. The Galerkin matrix elements B evolves according 
to 

B = -i[H,B]-{G,B}+2£r>c*Bct, ( 33 ) 



while the SPFs evolve according to 



klm 



(2) 
jklm< 



(34) 



where S jk, S^ lm and V] m {x) were defined in Eqs. |3()| >, |3lj 



and ( 32 1, respectively. 

Eq. \3A) is virtually identical to the standard ^-MCTDH 
equation of motion for the SPFs. The only difference lies in 
the definitions of S and S ^ . As for Eq. ( 33 I, we see that the 



main difference lies in the evolution of a matrix B instead of 
a coefficient vector. 

Eq. ( [34| is typically discretized using discrete-variable 
representation (DVR) techniques, FFT methods, finite dif- 
ferences or similar ifTTl . and the single-particle operator 
h — ir is then represented correspondingly. Equivalently, 
the single-particle space is approximated by the finite- 
dimensional space dictated by the discretization, inducing a 
finite-dimensional Fock space to begin with. Note that even 
though h — ir is non-Hermitian, orthonormality of tp is con- 
served during evolution. 

Like standard ^-MCTDH, the matrix S needs to be in- 
verted to evaluate the SPF differential equation; a well-known 
issue with MCTDH-type methods. It may happen that S be- 
comes singular for some reason, in which case a regulariza- 
tion approach is needed ifTTIl . In most applications, this hap- 
pens very rarely; typically at t = due to the choice of initial 
conditions, but experience suggests it does not affect the fi- 
nal results. This is, however, not trivial from a mathematical 
point of view, and for the sake of definiteness in the present 
work, we shall check that S is non-singular for our numerical 



experiment in Sec. IV 



Eq. p3) s hould be compared with the original Lindblad 
equation \12). Also, if tp = 0, we obtain the variational equa- 
tion of motion in a fixed linear basis, i.e., what is obtained us- 
ing a full configuration-interaction type approach. However, 



the Galerkin matrices defined in Eqs. ( 23 1 and ( 24 1 have time- 
dependent coefficients hm, ujkbn, and T jk which must be com- 
puted along the flow. This is a non-trivial task in general, and 
techniques common for ^P-MCTDH can be employed to deal 
with this in approximate ways lISl fTTIl . 

Note that B retains the natural block structure with respect 



to the number of particles, cf. Eq. ( 14 1. Eq. ( 33 1 can be written 

as 



B„ = -i[H,B„] - {G,B„} +2£r / *c*B B+1 ct 

and this is the most memory-economical representation, since 
the off-diagonal blocks vanish if p(0) is a pure state with N 



particles. In that case, B^ can furthermore be represented by a 
vector G % at all times, and there is no need to propagate 
the full block. Due to the presence of all the blocks in the 
definition of S and S* 2 ', however, this pure state cannot be 
evolved with ^P-MCTDH independently of the other blocks. 

If the dissipative terms vanish (i.e., if F — 0), and if p(0) 
is a pure state, the evolution is easily seen to be equivalent 
to a *P-MCTDH calculation. In p-MCTDH, tr(p) is not in 
general conserved 1171 . but it is so for closed systems, i.e., 
when £F(p) = 0. However, it is easily checked in the present 
case that, indeed, 

£tr(p) = Atr(B) = 0. 

Also, energy tr (Hp) is exactly conserved if tr(rp) = 0, that is 
to say, whenever the system does not touch the CAP. 

For the actual implementation of the evolution equations, 
it is useful to employ a generic enumeration scheme for the 
many-body basis states. In many-body codes, the Galerkin 
matrices (other than B) are rarely constructed in memory; 
instead the single- and double-particle integrals are kept in 
memory and the explicit action of H is computed using 
Eq. (23 i, for which the action of Cj and are implemented 
for example via mapping techniques as suggested in Ref. l30l : 
or simply using binary integers to represent a fermion state, 
and bitwise manipulations to define the action of c ; - etc.; a 
common technique in many-body nuclear physics calculations 

ED- 

As for choosing intitial conditions, we observe that as a 
generalization of ^-MCTDH capable of treating particle loss, 
a pure state p(0) = I'P/v) (^nI will be the usual choice. In that 
case, experience from ^P-MCTDH can be applied IfTTl . Typ- 
ical choices are single determinants/permanents or stationary 
states computed by imaginary time propagation, or combina- 
tions thereof as in the numerical experiment below. 



IV. NUMERICAL EXPERIMENT 

We present a numerical experiment for a model problem 
consisting of spin-polarized fermions in one spatial dimen- 
sion. We will study a situation where the initial state is a pure 
state with N = 3 particles whose norm gradually decreases 
due to a CAP. The situation is similar to the study in Ref. (9). 

We truncate the domain IR to [-R, +R], where R = 20. The 
single-particle Hamiltonian of our model is 



h = T + V(x) = -~ + V{x) 
where the one-body potential is of Gaussian shape 
V(x) = -86XPI-1.25X 2 ]. 



(35) 



Numerically we find that V(x) supports 4 bound one-body 
states. We choose a very simple CAP of standard power-form: 



r(x) = e(\x\-R')(\x\-R 



l\2 



(36) 
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FIG. 2: Single-particle trap potential V(x) and complex absorbing 
potential T(x) used in the numerical experiment. The trap poten- 
tial is Gaussian, see Eq. ^35\ , and the absorber is quadratic outside 
[-16,16], seeEq. {36}. 



where 8(x) is the Heaviside function. The particles are unaf- 
fected by the CAP in the region [—7?',/?'], where we we set 
R' = 16. We have verified that in the energy ranges of the cal- 
culations, very little reflection or transmission is generated by 
r. Fig.[2]shows the Gaussian well and the absorber. 

The particles interact via a smoothed Coulomb potential, 



u(xi,X2) — 2[(xj — X2) 2 + 0.1 2 



-1/2 



which is long-ranged. 

For discretizing the one-body space, we choose the stan- 
dard fast Fourier transform based method with N gr id = 128 
equidistant points with spacing Ax = 2/?/A^ gr id 11321 . 

For propagating the master equation, we choose a varia- 
tional splitting scheme [29 33 1, propagating the equations of 
motion with H' = T, i.e., kinetic energy only for a time step 
r/2, and then H' = H — T for a time step r, and finally H' = T 
for a time step r/2 again. This constitutes the propagation of a 
complete time step r. While being simple, the scheme has the 
advantage of having local error ^(t 3 ), that the r -propagation 
is numerically exact, and that the time step is not restricted to 
be t = ff(Ax 2 ). The potential step is integrated using a stan- 
dard explicit fourth-order Runge-Kutta method for simplicity, 
which is sufficient for our purposes. 

The initial condition is chosen as follows. Let x ¥2{x\,xi) 
be the two-body ground state of the CAP-less Hamiltonian. 
This is computed numerically by propagating the standard X V- 
MCTDH equations in imaginary time t = —is using L = 4 
single-particle states. It follows, that this state is also a sta- 
tionary for the present p-MCTDH method with a CAP as long 
as the overlap with the CAP is negligible. We have checked 
that this is indeed the case: Propagating the master equation 
with the two-body state as initial condition leads to an ab- 
sorption probability of 2.7 x 10~ 10 at t = ffinai = 30 which can 
safely be ignored. 

We act upon ¥2 with a creation operator c(g) J , where g(x) 
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FIG. 3: Space-time graph of square root ^/n{x,t) of particle den- 
sity in numerical experiment. Darker areas have higher density. The 
square root enhances contrast, but exaggerates low densities. From 
the plot, we can see that a single-particle function of Gaussian form 
is scattered off a bound two-particle state in a Gaussian well. The 
reflected and transmitted parts are absorbed by the CAP, revealing an 
oscillating trapped function of fewer particles. 



is a Gaussian of the form 

g{x) = eCexp[-(x + 2) 2 /0.75 + i3x], 

where Q projects away the 4 SPFs in the initial condition. g{x) 
describes an incoming particle of momentum k = 3 starting 
out at xq = —2. The final three-body initial state is then 

p(0)=J(g)\y 2 }(V 2 \c(g). 

The initial tp then consists of the L = 5 functions consisting of 
the 4 SPFs from the ground-state computation, and the single 
state g(x). 

Using this initial condition, we propagate p[t) for t < ffi na i = 
30. Fig. [3] shows a space-time graph of the particle density 
n[x,t) given by 

n(x,t) = tr[^{x)f{x)p{t)}. 

As expected, the plot shows the initial advance of the Gaus- 
sian wavepacket and its scattering off the well and the two- 
particle ground state. It is seen that scattering occurs both in 
the forward and backward direction. The scattered probabil- 
ity is absorbed upon entering the region |x| < R', and a system 
composed of less than three particles is seen to remain. Super- 
ficially, it is an oscillating two-particle system. The system's 
energy is E ss —7.355. 

However, the process is more complex, and by computing 
the probabilities p n (t) = tx\p n (t)] of having n particles in the 
system we may see what happens in more detail. In Fig. [4] 
p n is plotted for each 1 < n < 3. As the scattered probabil- 
ity density is absorbed, the probability of having n = 3 par- 
ticles decreases and the probability of n = 2 increases corre- 
spondingly. However, especially the absorption of the back- 
scattered wave reveals something interesting: the probability 
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FIG. 4: Plot of the probabilities p n = tr(p„) of having n particles in 
the system as function of t. The absorption events seen in Fig. [3] is 
seen to correspond to marked changes in the probabilities, p$ ls too 
small to merit an interesting plot. 



FIG. 5: Smallest eigenvalue of the matrix S which needs to be in- 
verted at each t. As one particle is almost entirely absorbed, the 
smallest eigenvalue falls off rapidly. It peters out at ~ 1.57 x 10~ 4 , 
safely away from zero. 



of having n = 1 particle in the system clearly becomes sig- 
nificant in this process: the bound two-particle system has a 
significant probability of being ionized by the collision, leav- 
ing a single particle. By inspecting the probability density 
ni(x,t) = tiffi (x)ifr(x)pi(t)] we verify that it corresponds to 
a bound one-body state superimposed on the two-body state. 
The probability po « 3.43 x 1CU 4 at t — fg na i, showing a very 
small probability of all particles vanishing. It is therefore not 
plotted. 

Although the initial bound two-particle state had negligible 
overlap with the CAP, there may still be errors introduced by 
a placing the CAP too close to the interacting system. For 
example, if a particle is absorbed prematurely, the remaining 
system may miss some correlations. Moreover, there seems 
to be a finite remaining probability of having three particles in 
the system. This is most likely due to reflections off or trans- 
missions through the non-ideal absorber F (which only was 
chosen for illustrative purposes), and not a bound three-body 
state. We have not investigated this in detail for the present 
experiment. 

Finally, we address the non-singularity of the matrix S in 
Eq. p4| ). In Fig. [5] the smallest eigenvalue of S is plotted 
as function of t. At the very last leg of the evolution, this 
eigenvalue cr m ; n drops off quickly. 

A sharp fall-off of cr m ; n is to be expected when a particle 
is almost entirely absorbed. This may be understood in terms 
of a non-interacting system. If interactions are not present, 
the whole system may be described by the SPFs alone, each 
evolving according to the non-Hermitian Schrodinger equa- 
tion under the single-particle Hamiltonian h — \Y. Thus, the 
eigenvalues cr; only change because of the CAP. As a particle 
is absorbed, one eigenvalue goes to zero. 

Usually, small eigenvalues crj give very rapidly changing 
SPFs. However, we have observed that the natural orbital cor- 
responding to cr min , which we observe resides in the CAP 
region and thus represents the absorbed particle, does not 



change significantly in the last leg of the evolution. This indi- 
cates that the SPF no longer becomes relevant for the descrip- 



tion in in the sense that the right-hand side of ( 34 1 de-couples 
from this SPF. In this way, the near-singularity of S stemming 
from particle absorption may in fact have no impact at all on 
the evolution. This may be related to the simple fact that for a 
pure two-fermion system, the single-particle reduced density 
matrix is singular whenever L is an odd number. (It is easy to 
show that the eigenvalues are zero or of multiplicity 2 in the 
two-fermion case.) We cannot draw any firm conclusions con- 
cerning this from our simple experiment, except for pointing 
out that S becoming singular may have different causes and 
consequences compared to pure state-MCTDH. 



V. MIXTURES OF SPECIES 

The density operator MCTDH formulation for identical 
fermions is readily generalized to mixtures of arbitrary num- 
ber of species of particles, such as mixtures of 3 He and 
4 He (fermions and bosons). Using second quantization, the 



derivation becomes analogous to the treatment in Sec. Ill and 
Ref. l27l . so we only state the main results here. 

As Ref. E71 we consider two different species of particles A 
and B for simplicity, as the generalization to K species follows 
immediately. Each species have a Fock space Jj?® , i — A,B. 
The total Hilbert space is the product space 

Each species is assigned a set ip( l > of single-particle states 
ip^(x) and ip^ B \y), but they have no a priori connection, 
as the single-particle spaces may be very different. Conse- 
quently, the operators a.j (for species A) and typ (for species 
B) all commute since the species are distinguishable from 
each other. As previously, the creation operators are used to 
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construct finite-dimensional Fock spaces T^W with determi- 
nant or permanent basis functions <5>^ . For the product space, 
the basis functions are 



(A), 
J[n] 



Note that as the particles are distinguishable, we speak of 
(n,m)-particle states. Fock space is divided into subspaces 
with n particles of species A, and m particles of species B. The 
density operator p will then be block diagonal with respect to 
the particle numbers: 

N A N b 
n=0m=0 

where Nj are the maximum number of particles in the system, 
determined by the initial condition. 

Each species has its internal Hamiltonian, but for the equa- 
tions not to separate into the previously studied case, we need 
an interaction. A generic two-body inter-species interaction 
may be written as 

n m LW L<- b > 

W = £ Y, w ( x i>yj) = E E w jkim.a)aib\b m 

i= 1 j= 1 /,/= 1 k,m— 1 

in first and second quantization form, respectively. Here, 

/ (A) (B) I / x, (A) (B)\ 

The usual factor 1 /2 is not present, since the particles are not 
identical. 



Each species also has its own absorber rW, which need not 
have any a priori relation. 

Working through the equations of motion, noting that each 
species' SPFs are independent from each other, we obtain the 
following equation for the Galerkin matrix blocks B„ ,„: 



B„ 



= -i[H,B„, m ]-{G^+G^ B„, m } 

- 2 L r X ) ^ B «+i.'» a /+ 2 L r ? b 't B «.'»+i b ; t 

jk jk 



with H = H* A ) +H' b ' +W and an otherwise obvious notation. 

We obtain an SPF equation of motion for each species. 
They will contain species-specific analogues of S, and 
of the mean fields C7® of U, and also inter species analogues 
of the reduced two-body density matrix elements S ( AB - 2 ) and 
mean fields due to W, exactly as in the 'P-MCTDH |27l . 
These are defined by 



, (AB,2) _ c (BA,2) _ 



jklm 



kjml 



ti{p 2 a]aiblb m ), 



and 



Wg> = / <P^(x)w(x,y)^ A \x)6x, 



respectively. We get the SPF equations of motion 



L. (A) „ (A) 

j 

L. (B) „ (B) 
Vj S k/ 



kj 



Q {A) 

Q(B) 



j klm klm 



klm 



klm 



These equations have an obvious symmetry with respect to 
particle species interchange. These equations are identical to 
those given in Ref. Il27l . except for the reduced density matri- 
ces being defined in terms of a density operator and not a pure 
state. 

The generalization to K species is straightforward, and we 
refrain from going into further detail. Note however, that 
one interesting special case is obtained when the number of 
species equals the number of initial particles, so that the ini- 
tial state is a density operator in the space 

JT = M^ l) ® ^ (2) ® • • • ® J^I K) . 

In this case it seen that, as all particles in the system are in fact 
distinguishable, we have obtained the usual MCTDH method 
for density operators, but with a CAP. Of course, the same can 



be said of pure-state MCTDH for mixtures of particles - as we 
approach K = N species, where is the number of particles, 
we are back at plain MCTDH for distinguishable particles, 
and the circle is closed: MCTDH for identical particles can be 
viewed as MCTDH with (anti-)symmetry constraints on the 
coefficients, and plain MCTDH can be viewed as A^-species 
MCTDH for mixtures. 



VI. CONCLUSION 

A system of jV particles described by a Hamiltonian with 
a complex absorbing potential evolves irreversibly in time. 
In order to describe the remaining particles as some are lost 
to the absorber, a master equation on Lindblad form in Fock 
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space is needed, as first demonstrated in Ref. \22\. This equa- 
tion was discussed at length, and a multiconfigurational time- 
dependent Hartree method was presented that is a strict gener- 
alization of standard pure-state MCTDH evolution for identi- 
cal particles or mixtures. A numerical experiment on a simple 
system of N = 3 spin-polarized fermions was reported. 
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